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ABSTRACT 

Wc  develop  a  semi-discrete  approximation  framework  for  linear  non- 
autonomous  nonhomogeneous  functional  differential  equations  of  retarded  type. 
The  approximation  schemes  are  constructed  and  convergence  results  are  obtained 
through  the  approximation  of  an  associated  abstract  evolution  operator. 
Computer  implementation  of  the  schemes  is  outlined  and  a  spline  based 
method  included  in  the  framework  is  constructed.  Extensions  of  the  semi¬ 


discrete  methods  to  schemes  incorporating  full  discretization  and  difference 
equation  approximation  are  also  discussed.  Numerical  results  for  several 
examples  demonstrating  the  feasibility  of  the  schemes  are  presented. 


Part  of  this  research  was  carried  out  while  both  authors  were  Visiting  Scientists 
at  the  Institute  for  Computer  Applications  in  Science  and  Engineering,  NASA  Langley 
Research  Center,  Hampton,  VA  under  NASA  Contracts  No.  NAS1-15810  and  NAS1-16394. 

^Work  supported  in  part  by  the  Air  Force  Office  of  Scientific  Research  under 
Contract  AFOSR  76-3092D,  in  part  by  the  National  Science  Foundation  under 
Grant  NSF-MCS  7905774-02,  and  in  part  by  the  U.  S.  Army  Research  Office  under 
Contract  ARO-DAAG29-79-C-0161. 


1. 


Introduction 


Our  presentation  deals  with  spline  based  approximation  theory  for 
linear  nonautonomous  nonhomogeneous  functional  differential  equations  (FDE) 
via  approximation  of  the  associated  evolution  operator  for  the  homogeneous 
equation.  Our  theoretical  framework  and  the  resulting  convergence  results 
are  analogous  to  the  Trotter-Kato  type  approach  for  convergence  of  solution 
semigroups  for  autonomous  delay  systems  as  developed  in  [8].  We  restrict 
our  discussions  to  piecewise  linear  spline  approximations,  but  have  attempted 
to  state  the  theoretical  results  in  a  form  so  that  one  might,  if  one  so  de¬ 
sired,  easily  extend  the  ideas  in  a  rather  straight  forward  manner  to  treat 
higher  order  spline  approximation  schemes. 

One  might  also  view  the  discussions  here  as  one  concrete  realization 
of  the  evolution  operator  theory  presented  for  abstract  nonautonomous  equa¬ 
tions  in  [14].  Indeed,  as  we  shall  see  below  our  approach  is  very  closely 
related  to  the  efforts  of  Reber  in  [17,  18]  who  used  the  Krein  approach  to 
develop  for  nonautonomous  FDE  control  problems  an  approximation  theory  based 
on  the  so-called  "averaging"  approximations  of  [4,5].  However,  since  the 
spline  based  methods  of  [8]  have  been  shown  to  offer  considerable  improve¬ 
ment  over  the  averaging  schemes,  and  since  the  needs  for  schemes  to  treat 
nonautonomous  problems  are  rather  obvious  in  a  number  of  areas  of  applica¬ 
tions  (for  one  example,  see  the  discussion  of  the  tracking  models  in  [24]), 
we  feel  that  our  modest  contribution  below  towards  the  development  of  spline 
schemes  is  warranted.  We  demonstrate  the  computational  efficacy  of  our 
ideas  by  presenting  a  sample  of  our  numerical  findings  in  section  6. 

Other  authors  have  considered  approximation  schemes  for  nonautonomous 
delay  systems.  Both  Delfour  [11]  and  Reber  [18]  consider  optimal  control 
problems  for  linear  nonautonomous  retarded  FDE  and  develop  "full  discretization" 
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techniques  (state  and  time  discretization)  that  employ  the  averaging  ideas 
of  [4,5] .  Kappel  and  Schappacher  in  [13]  consider  nonlinear  nonautonomous 
delay  systems  which  they  approximate  by  linear  interpolating  spline  schemes 
in  the  "state"  space  C.  This  results  in  essentially  an  averaging  type 
scheme  as  opposed  to  the  spline  schemes  discussed  below.  In  a  similar  spirit 
Kunisch  [15]  discusses  the  averaging  and  essentially  equivalent  linear  inter¬ 
polating  spline  schemes  for  optimal  control  of  neutral  FDE  with  nonautonomous 
retarded  terms. 

It  is  possible  to  take  other  approaches  to  a  spline  approximation  theory 
for  delay  systems.  In  [2]  we  combine  dissipativeness  of  nonlinear  operators 
with  Gronwall  estimates  to  develop  one  approach  to  approximation  of  nonlinear 
FDE.  These  ideas  are  pursued  in  [7,  10]  to  yield  spline  approximation  methods 
(in  the  context  of  parameter  estimation  schemes)  for  quite  general  nonlinear 
nonautonomous  systems  of  FDE.  In  addtition  to  differing  in  spirit  from  the 
approach  of  [7,  10],  our  results  below  are  applicable  to  linear  systems  with 
nonsmooth  coefficients  while  the  theory  of  [7,  10]  requires  some  smoothness 
on  the  coefficient  matrices  if  it  is  applied  to  linear  systems. 

For  our  detailed  development  below  we  consider  ordinary  differential 
equation  (ODE)  or  semi-discrete  approximations  to  linear  FDE.  These  ODE 
for  the  "Fourier"  coefficients  of  the  approximate  solution  relative  to  a 
fixed  spline  basis  must  then  be  solved  by  a  high  order  ODE  solver  (a  fourth 
order  Runge-Kutta  in  the  case  of  the  examples  presented  in  section  6).  An 
alternative  approach  invol\ ing  Immediate  full  discretization  (in  both  time 
and  state)  and  resulting  in  a  difference  equation  approximation  of  the  FDE 
could  be  taken  in  the  spirit  of  the  efforts  of  Rosen  in  [19,  20].  We  give 
a  brief  indication  of  some  of  these  ideas  in  section  5  but  will  not  pursue 
a  full  detailed  development  along  these  lines. 


— _ _ _ _ __  .  J 
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While  we  shall  not  discuss  such  applications  directly,  the  reader  should 
be  aware  that  the  approximation  ideas  developed  here  are  readily  (and  profit¬ 
ably)  used  in  optimal  control  and  parameter  estimation  problems  (e.g,  see  [1], 
[3],  [5],  [6],  [91,  [12],  [16],  and  [21]). 

The  presentation  below  is  organized  in  the  following  manner.  In  section 
2  we  summarize  equivalence  results  between  FDE  and  abstract  evolution  equations. 
VJe  also  establish  a  dissipative  condition  for  the  operators  involved  that  is 
crucial  in  our  development.  Basic  approximation  results  are  given  in  section 
3,  first  for  systems  with  continuously  differentiable  coefficients  and  then 
for  systems  with  coefficient  matrices.  This  is  followed  by  section  A  in 

which  we  develop  a  particular  piecewise  linear  spline  scheme  in  detail  and 
explain  how  it  is  to  be  implemented.  A  brief  discussion  of  full  discretization 
ideas  is  given  in  section  5.  We  conclude  the  paper  with  some  representative 
numerical  findings  in  section  6. 

Notation  throughout  is  completely  standard  with  respect  to  symbols  for 

C,  Lp,  etc.  We  denote  the  usual  Sobolev  spaces  W^^  of  functions  f  with 

(k-1)  (k)  k 

f  absolutely  continuous  and  f  in  Lj  by  H  .  For  Lebesgue  spaces 

of  Rn-valued  functions  on  (a,b)  we  adopt  the  notation  Ln(a,b)  while  L 

p  nxn 

denotes  the  space  of  n  square  matrices.  Finally,  we  shall  sometimes  use 
D<f>  to  represent  the  derivative  of  a  function  <f>. 


2.  The  Linear  Non-autonomous  Functional  Differential  Equation  and  its 
Equivalent  Formulation  as  an  Abstract  Evolution  Equation 

In  this  section  we  describe  the  functional  differential  equation  (FDE) 

initial  value  problem  for  which  we  seek  approximation  schemes  and  give  an 

equivalent  formulation  of  it  as  an  abstract  evolution  equation  in  an  infinite 

dimensional  Hilbert  space.  Many  of  the  results  and  ideas  which  are  outlined 

and  summarized  below  are  discussed  in  detail  in  [18,  sections  2  and  3]. 
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We  consider  n-vector  systems  of  the  form 

(2.1)  x(t)  =  L(t)xt  +  f(t)  t  >  0  , 

(2.2)  x(s)  =  n  xg  =  <j>, 

where  f  €  l”  ^oc(s,+°°),  ^  ^  rII»  #  e  L2(-r,0)  and  x£  denotes  the 
function  0  -*•  x(t  +  0),  -r  _<  0  <_  0.  We  assume  that  for  each  t  0,  the 
linear  operator  L(t) :L2(-r,0)  -*  Rn  has  the  form 

v  0 

(2.3)  L(t)<j>  =  £  A,(t)<K-T.)  +  f  A(t,e)4>(0)de, 

i=0  1  -r 

with  0  =  x0  <  <  t2  •••  <  Tv  =  r,  A1  €  ^([s,-^)  ,1^) ,  i  =  0,1,2  •••  V 

and  the  function  t  -*■  A(t,*)  an  element  of  L  ([s,-H»),  L~([-r,0],L  )).  The 

00  z  n*n 

point  evaluations  of  <{>  €1  L2(-r,0)  required  in  the  evaluation  of  L(t)<p 
pose  no  essential  conceptual  difficulties  since,  roughly  speaking,  we  shall 
interpret  solutions  of  (2.1)  as  functions  satisfying  that  equation  in  inte¬ 
grated  form  (i.e.  the  differentiated  form  in  the  almost  everywhere  sense)  and 
thus  any  occurance  of  L(t)<J)  for  <j)  only  an  L2  "function"  can  be  considered 
as  appearing  under  an  integral  with  the  L(t)(J>  then  denoting  an  equivalence 
class  of  functions.  For  a  further  discussion  of  this  point,  we  refer  the 
reader  to  [8]. 

A  solution  to  (2.1),  (2.2)  is  a  function  x:  [s-r,T]  -*■  Rn,  T  >  0,  such  that 
x  €  H^(s,T),  x  satisfies  equation  (2.1)  almost  everywhere  in  [s,T],  x(s)  =  n 
and  xg  =  (J).  Standard  arguments  [17]  can  be  used  to  show  that  the  FDE  initial 
value  problem  (2.1),  (2.2)  has  a  unique  solution  which  depends  continuously 
upon  the  initial  conditions  and  the  nonhomogeneous  term  f.  We  shall  on 
occasion  employ  the  notation  x(t;i|,<j>,f)  (and  (rj ,<ji, f ) )  in  order  to  denote 
this  unique  solution  (and  its  past  history  on  [t-r,t])  to  (2.1),  (2.2)  corre¬ 
sponding  to  a  particular  choice  of  n»<f>,  and  f. 


-5- 


We  next  define  the  Hilbert  space  Z  by 

Z  =  Rn  *  Lj (_r ,0) , 

with  inner  product 


and  reformulate  the  FDE  initial  value  problem  (2.1),  (2.2)  as  an  abstract 
evolution  equation  in  Z.  Corresponding  to  f  =  0,  it  is  possible  to  define 
a  solution  operator  for  (2.1),  (2.2)  on  Z  by 


U(t,s)  (n,<fi)  =  (x(t;n.<K0)  ,xt  (n,4>,0)) . 


It  is  easily  verified  that  for  T  >  s  the  operators  (U(t,s):t  €  [s,T]} 
are  continuous  in  t  and  uniformly  bounded.  In  addition  U(t,s)  is  an  evolution 
operator  on  Z  in  that  the  uniqueness  of  solutions  to  (2.1),  (2.2)  guarantees 
that  it  satisfies  U(s,s)  *  I  and  the  transition  property  U(t,s)  =  U(t,T)u(T,s) 
for  all  s  <  T  <  t  <  I. 

Returning  to  the  nonhomogeneous  problem,  we  define  for  arbitrary 
f  €L2,1oc<8’4~) 

t 

(2.4)  z(t,s;n,4>,f)  =  U(tjs)  (n,<J>)  +  J  U(t,s)(f  (s),0)ds. 

S 

For  each  (n,4>)  C  Z,  f  €  lo^3*"1"00)  and  t  s  the  expression  given  in 

(2.4)  exists  and  is  continuous  in  t.  Furthermore,  it  can  be  shown  that 


(2.5)  z(t,s;n,<f>,f)  =  (x(t;n,ij>,f),xt(n,<J>,f)). 

Equation  (2.5)  states  that  (2.4)  and  (2.1),  (2.2)  are  equivalent,  and  in 
fact,  forms  the  basis  for  the  approximation  schemes  developed  in  the  next 
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section.  Indeed,  we  construct  convergent  approximations  to  the  solu¬ 
tion  of  (2.1),  (2.2)  via  the  construction  of  convergent  approximations 
to  z  of  (2.4). 

If  one  is  willing  to  impose  additional  restrictions  on  (2.1),  (2.2)  a 
stronger  result  can  be  established.  Consider  the  initial  value  problem 
(2.1),  (2.2)  with  the  alditional  assumption  that  the  coefficient  matrices, 
the  kernel  in  the  distributed  term  and  the  nonhomogeneous  term  be  continuously 
differentiable  in  t  and  that  (t|><J>)  €  W  =  {(r|,<j>)  6  Z:<J>  €  H^(-r,0)  ,n  =  0(0)  }. 
It  then  can  be  shown  that  z  given  in  (2.4)  is  the  unique  solution  to  the 
abstract  evolution  equation  in  Z  given  by 

(2.6)  z(t)  =  A(t)z(t)  +  (f(t),0), 

(2.7)  z(s)  =  (0,40  » 

where  for  each  t  >  s  the  operators  A(t)  :W  CZ+Z  are  defined  by 

(2.8)  A(t)(4>(0),4>)  =  (L(t)<}>,4) . 

In  addition  it  can  be  verified  that 

w(t)  =  (x(t;n,4>,f) ,  xt(o,4>,f)), 

is  also  a  solution  to  (2.6),  (2.7)  and  hence  must  coincide  with  z.  Thus 
under  these  stronger  hypotheses  (2.1)  -  (2.2),  (2.4) ,  and  (2.6)-  (2.7)  are  all 
equivalent . 

The  existence  of  an  inner  product  on  Z,  equivalent  to  the  standard 
inner  product  on  Z  defined  above,  and  an  w  for  which  the  operator 
A(t)  -  <-'I  is  dissipative  plays  an  essential  role  in  many  of  the  arguments 
which  follow.  Toward  this  end,  we  define  the  same  inner  product  on  Z  as 
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that  one  employed  In  [8]  and  [10]  for  similar  purposes.  Let  p  be  the 
step  function  on  [-r,0]  defined  by 

P(6)  =  j  "Vj+1  -  6  <  "Vj*  j  =  v’ 

and  Z.(  the  space  Z  with  inner  product  <•,♦>  are  given  by 

t  0 

<(£,<Jj),(n,<f>)>  =  Cn  +  f  <i>(9)<t>(0)p(0)dd  . 

P  -r 

It  is  easily  verified  that  the  <•,•>  inner  product  is  equivalent  to 
the  standard  inner  product  on  Z  and  moreover,  the  following  lemma  can  be 
established. 

LEMMA  2.1.  For  each  t  >  s,  A(t)  -  wl  is  dissipative  in  Z  .  That  is 

P 

<A(t)z,2>p  _<  W<Z,Z>p, 

for  each  z  £  W  with 

3  .v  2  .  V  .  , 

W=2ml  +  2ml+2  +  1 

and 

v 

®i  =  E  iAilco +  lAL  • 

i=0 

It  is  in  fact  the  case  that  ^(t)  -  wl  is  maximal  dissipative.  That  is  to 
say,  A(t)  -  w I  is  onto.  The  verification  of  this  latter  claim  can  be  found 
in  [18,  section  3].  The  reader  should  note  that  while  the  hypothesis  of 
Lemma  3.3  in  [18]  include  the  assumption  of  smooth  coefficients  on  the  right 
hand  side  of  the  FDE,  it  is  easily  seen  that  this  assumption  does  not  play  a 
role  in  the  arguments  used  to  show  that  A(t)  -  wl  is  onto. 

RemcUik  2.1.  By  the  Lumer-Philllps  theorem  (cf.  [25,  section  IX. 8]) 
the  fact  that  for  t  €  (s,+°°)  fixed,  A(t)  -  wl  is  a  maximal  dissipative 
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operator  on  the  Hilbert  space  Z  is  sufficient  to  conclude  that  A(t)  is 
the  infinitesimal  generator  of  a  Cq  semigroup  of  bounded  linear  operators 
on  Z. 

3.  Approximation  Results 

Following  the  ideas  discussed  in  [8]  we  base  our  approximation  schemes 

on  the  following  construction.  For  each  N  =  1,2,  •  •  • ,  {Z.,,P,,,A  (t)  }  will  be 

N  N  N 

called  an  approximation  scheme  for  the  initial  value  problem  (2.1),  (2.2)  if 
{ Z^j }  is  a  sequence  of  finite  dimensional  subspaces  of  Z^,  { }  is  a  sequence 
of  operators,  where  for  each  N,  P^:Zp  Z^  is  the  orthogonal  projection  from 
Zp  onto  Z^  and  (A^Ct)}  is  a  sequence  of  t-dependent  operators  on  Z^. 

THEOREM  3.1.  Consider  the  FDE  initial  value  problem  (2.1),  (2.2)  under 

the  additional  assumption  that  A.  £€'*'(  [s ,+°°)  ,  L  )  ,  t  -*■  A(t,«)  £C^([s,+°°), 

i  nxn 

^2 ( (-r,0) , LnXn) )  and  suppose  that  {Z^,P^,A^(t) }  is  an  approximation  scheme 
for  (2.1),  (2.2)  satisfying  the  following  conditions 

(1)  Z^  CW  =  dom(A(t))  N  =  1,2,***  . 

(2)  For  each  t  >^s,  A^(t):ZN  ZN  is  defined  by  AN(t)  =  PNA(t),  N  =  1,2 

(3a)  lim  P  z  =  z  in  Z  for  all  z  £  Z. 

N-*°° 

(3b)  For  ip  €  W  with  P^ip  =  PN(ip(0),ip)  =  (ipN(0),ipN)  we  have 

lim  L(t)ip  =  L(t)tp  in  Rn  for  each  t  £  [s,T]  and 
N~x»  ^ 

lim  Dtp.,  =  Dv(j  in  L”(-r,0)  with  |D(tp.,  -  tp)  j  <  K|D2tp|,  K  indepen- 

N-KX)  n  1  N 

dent  of  N  and  ip  for  all  tp  €  W  with  tp  £H2(-r,0). 

Then  if  UN(t,s)  denotes  the  evolution  operator  (fundamental  matrix  solution) 
corresponding  to  the  finite  dimensional  ordinary  differential  equation  in  Z^ 


given  by 
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z^(t)  -  A^(t)z^(t), 

we  have  that 


li»  |[PNU(t,S)-Vt,S)P[IJz0|  -  0, 
for  each  zQ  G  Z,  uniformly  in  t  for  t  G  [s,T]. 


Vfioof,.  An  application  of  Lemma  2.1  and  the  fact  that  the  P„  are 
-  N 

orthogonal  projections  yields  the  following:  For  €  Z^ 

</^N(t)zN’ZN>p  =  <PNA(t)zN,ZN>p  =  'A(t)zN’PNZN''p 


=  <A(t)zN,zN>p  <  «  <zN>  VP> 


where  u)  is  defined  in  the  statement  of  Lemma  2.1,  and  is  independent  of  N. 
The  calculation  above,  and  the  fact  that  /^(t)  is  a  bounded  linear  operator 
defined  on  the  finite  dimensional  space  Z^  are  sufficient  to  guarantee  that 
Y*>  -  ojI,  N  =  1,2,***,  are  maximal  dissipative,  and  moreover 

cKA^t))  C  {AC<{::ReA  <  oj},  N  =  1,2,**«,  t  >  s. 

Thus  for  all  A  6  ij;  with  Re  A  w  the  resolvent  operators  R(A,A^(t))  exist 
and  by  standard  arguments  (c.f.  [14,  p.  85])  we  have 

(3.1)  |R(A,yt»|p  =  KYO-AI)*1^  1  N  =  1,2,  •  •  • ,  t  €  [s ,+°°)  . 

Inequality  (3.1)  and  the  same  arguments  used  to  establish  the  validity  of 
Theorem  3.5  in  [18]  allows  us  to  conclude 

|UN(t,s)|z  <  Mew(t“s), 

where  M  is  a  constant  independent  of  N  and  L(*)>  the  homogeneous  part  of 
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the  right  hand  side  of  the  FDE  (2.1).  The  reader  is  instructed  to  note 
that  while  the  constant  M  derived  in  the  proof  of  Theorem  3.5  in  [18] 
does  depend  on  L('),  our  M  does  not.  This  is  a  consequence  of  the 
fact  that  the  weighting  function  p  in  the  inner  product  <•,•>  *-s  ^n“ 

dependent  of  the  A.,  whereas  that  in  the  inner  product  <♦,•>.  chosen  by 
Reber  is  not.  The  significance  of  our  choice  for  the  inner  product  will 
be  apparent  in  the  discussions  below  pertaining  to  extensions  of  our  con¬ 
vergence  results  to  equations  with  the  A^,  A  only  in  L^. 

Let  D  be  the  subset  of  Z  defined  by 

d  =  { (n,4>)  ez:4>  ec2(-r,0),n  =  <K0),L(s)<j>  =  $(o)}. 

Using  the  fact  that  A(s)  is  the  infinitesimal  generator  of  a  semi¬ 

group  of  operators  on  Z  (cf.  Remark  2.1)  arguments  similar  to  those 
used  to  verify  Lemma  2.2  of  [8]  can  be  used  to  establish  the  fact  that  D 
is  dense  in  Z.  Consider  next  the  initial  value  problem  in  Z  given  by 

(3.2)  2 ( t)  =  A(t)z(t), 

(3.3)  z(s)  =  zQ  ,  zQ  €  D, 

and  the  following  identity  derived  from  it: 

V(t)  =  Vt)PNz(t)  +  tPNA(t)  Z(t>  ~  AN<t>pNz(t)  i - 

Recalling  that  the  are  orthogonal  projections,  and  therefore  uniformly 

bounded,  we  have 

£  [PNz(t)]  =  AN(t)[PN2(t)]  +  [PNA(t)2(t)-AN(t)PNz(t)], 


and  thus,  by  the  variation  of  constants  formula 
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(3.4)  P  z(t)  =  UN(t,s)PNz(s)  +  f  UN(t,T)[PNA(T)  -  AN(T)PN]z(T)dT. 

s 

Since  Zq  G  DCW,  the  unique  solution  z  to  (3.2),  (3.3)  is  given  by 

z(t)  =  U(t,s)zQ, 

and  hence  (3.4)  can  be  rewritten  as 

PkU(t,s)z0  -  Vt’s)PNZ0  =  /UN^»T)fV(T)"AN(T)PN]z(T)tiT- 

s 

Therefore 

(3.5)  |P  U(t.s)z0-U  (t,s)PNz  I  -  1/  UN(t,T)[PNA(T)-AN(T)PN]z(T)dxr 

s 

1  /  IVt,T>  I  I  [V(T)  "  \(t)Pn]z(t)  IdT 

s 

<  Mew(T's)  f  |fPNA(T)-AN(T)PNlz(T)|dT. 
s 

2 

Let  z  »  (<p(0),<p)  CZ  be  such  that  <p  e  H  (-r,0). 

Then 

|  (PNA(t)  -  AN(t)PN]z  1  2  =  |[PNA(t)-PNA(t)PN](4>(0),1l))|2 

<  |  [A(t)  -  A(t)PN3  (0(0)  ,4>)  | 2 

=  |  (L(t)$,D<j>)  -  (L(t)ct»N,D4)N)  I  2 

=  |(L(tH-L(t)<|>N)|2n  +  |  D<b  -  D<bN  i  J  -0, 

R  2 
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as  N  -*•  »  for  each  t  6  [s,T]  by  condition  (3b).  Thus  P.  A(t)z  -  A„,(t)P.,z  >  0 

N  N  N 

for  each  t  G  [s,T],  and  each  z  =  C<t (0)  »4>)  G  Z  with  <f>  GH2(-r,0). 

Now  under  the  smoothness  assumption  on  L(*)  for  zQ  =  (4>(0 ),<}>)  G  D 

2  9 

we  have  that  x(- ;  <f>  (0)  ,<M)  6H  (s,+“>).  Consequently  x  (<J>(0)  6H  (-r,0) 

for  t  >_  s.  Therefore  in  the  light  of  the  equivalence  established  in  section 

2,  and  the  arguments  above,  we  have 

(3.6)  |[P/(t)  -AN(t)PNlz(t)|  =  |[PNA(t)-AN(t)PNJU(t,s)z0| 

=  |  tP„A(t)  -  AH(t)PN](x(t;0(O),«,O),xt(.KO),«ji,O))| 

0  for  each  tS[s,T]  with  the  convergence 

2 

being  dominated.  Indeed,  x  c  H  (s,+")  and  an  application  of  condition  (3b)  yielda 

|  [P/(t)  -  AN(t)PN]z(t)  |2  <  |  [A(t)  -  A(t)PN]  (x(t;<J>(0)  ,<p,0)  ,xt(<K0)  ,4>*0))  I  2 

=  i(L(t)xt(<K0),<h0),Dxt(<K0),<l>,0)) 

-  (L(t)x^(<}>(0) » cp ,  0)  ,Dx^(<})(0) , , 0) )  |  2 

=  |  L(t)(x  (<K0M,0) -x%(0),<N0))|2 

R 

+  |D(xt(*(0),*,0) -x^(<t>(0),<t,0))|2 

<  |L(t)(x  (4»(0)  ,4>,0) -x^(c(>(0),lj),0))|2  +  K2  |  D2x  (*«», 

R  i 

J 

<  |L(t)(x  ((K0),*,0)  -  x^C0(O)  ,d>,0))  )  2n  +  K2|D2x|2n  | 

R  L^s-r, 

where  the  bound  on  the  second  term  in  the  above  estimates  is  a  consequence  of 

N  N 

condition  (3b)  and  (x  (t),xt)  =  P^(x(t),xfc).  We  note  that  it  is  not  true  in  I 
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N  N 

general  that  x  ■  (x  )fc.  To  see  that  the  convergence  in  the  first  term  is 

N 

dominated  as  well,  we  first  observe  that  x  (t)  -*•  x(t)  uniformly  in  t  for 
t  £  [s-r,T].  This  follows  from  condition  (3a)  and  the  fact  that  the  subset 
of  Z  given  by  (z(t):t  e.  [s,T]}  is  compact  (being  the  continuous  image  of 
a  compact  set  in  R) .  Then,  since  for  -r  £  0  £  0 


we  have 


x!?(0)  =  x^(0)  +  f  (Dx^) (a) da, 
t  t  jQ  t 


|xj(0)-xt(0)|  £  |x^(0)  -  xt(0)  |  +  /|D(xJ-xt)(a)|da 

0 

i  |xj(0)-xt(0)|  +  r^|D(x^-xt)|2 


£  |xN(t)  -x(t)|  +  A|D2xt|2 


£  |xN(t)-x(t)|  +  r^K|D2x 

L2(s-r,T) 


£  |xN-xl  +  A|D2x 


L2(s-r,T) 


Therefore 


|  L(t)  (xt  (<f>(0)  ,cf>,0)  -  x^(0(O)  ,0,0))  |  n 

R 

<  (|xN(*;0(O),0,O)  -  x(* ;0(O) ,0,0) |  +  r^K|D2x|  \ 

'  Lj(s-r,T) > 

•  (Ely*>L  +^!(  /|A(-,0)|2de)4|a>|. 


Returning  to  (3.5),  the  convergence  stated  in  (3.6)  together  with  the  dominated 

convergence  theorem  allow  us  to  conclude  that  for  each  Zq  £  D 

PNU(t,s)zQ  -  UN(t,s)PNZp  |  +  0  as  N  -*■“>,  uniformly  in  t  for  t  6  [s,T].  But 


D  is  a  dense  subset  of  Z  and  the  operators 
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[PNU(t,s)-UN(t,s)PN]  , 

are  uniformly  bounded  in  N.  Therefore 

|PNU(t,s)z  -  UN(t,s)PNz|  ->0,  N  -*■  °°, 

uniformly  in  t  for  t  6  [s,T]  for  each  z  6  Z. 

RejmaAk.  It  is  possible  to  obtain  estimates  for  the  rate  of  convergence 
of  the  evolution  operators  on  restricted  classes  of  initial  data.  The  rele¬ 
vant  arguments  are  in  the  same  vein  as  those  used  to  derive  the  estimates 
for  the  rate  of  convergence  in  the  autonomous  case.  The  details  of  these 
arguments  can  be  found  in  [8]. 

Turning  our  attention  to  the  nonhomogeneous  problem,  for  (o>4)  6  Z 
and  f  €  l”  we  consider  the  approximating  finite  dimensional  non¬ 

homogeneous  ordinary  differential  equation  in  Z^  given  by 

Vfc)  =  Vt)zN(t)  +  pN<f(t)*0)» 

ZN(S)  =  V1’40* 

Using  the  variation  of  constants  formula,  we  can  write  down  the  solution 

to  this  initial  value  problem.  It  is  given  by  zN(t,s;n»4>f)  =  UN(t,s)PN(o,4) 
t 

+  J U^(t ,o)PN(f (a) ,0)da  for  t  ^  s.  An  application  of  arguments  analogous  to 
s 

those  used  to  verify  Theorem  3.2  of  [8]  will  establish  the  validity  of  the 
following  theorem. 

THEOREM  3.2.  Under  the  hypotheses  and  conditions  of  Theorem  3.1  we  have 

(a)  For  (0,4))  6  Z  and  f  e  L^CsjT) 

lim  z  (t,s;n,4»f)  =  z(t,s;u,4>f) , 

N-**> 


uniformly  in  t  for  t  €  [s,T]  and  uniformly  in  f 
for  f  in  bounded  subsets  of  L^s.T). 

(b)  For  (x^O.y^t))  6  ZN  defined  by  (x^t)  ,yN(t))  =  zN(t,s;  n,  <M)  , 
we  have 

lim  x^(t)  =  x(t;n,<))»f), 

N-+® 

uniformly  in  t  for  t  €[s,T]. 

(c)  If  { f ^ }  is  a  sequence  in  L^Cs.T)  converging  weakly  to  f  then 

lim  zN(t,s;n,<J),fk)  =  z(t,s;n,<M). 

N,k-*>° 

uniformly  in  t  for  t  €[s,T]. 

RejnaAk.  Although  it  will  not  be  discussed  in  this  paper,  part  (c)  of 
Theorem  3.2  above  plays  an  essential  role  in  the  application  of  our  approxi¬ 
mation  results  to  the  solution  of  optimal  control  problems  governed  by  FDE 
of  the  form  (2.1).  These  ideas  are  discussed  in  detail  for  the  case  of  an 
autonomous  equation  in  [5]. 

We  conclude  this  section  with  a  discussion  of  the  details  involved  in 
extending  the  approximation  results  above  to  FDE  initial  value  problems  of 
the  form  (2.1),  (2.2)  with  non-smooth  right  hand  sides.  Define 

A  -  [.,,(<=,+.),(  X  !„„„)*  1 -2(<-r.°>  •!„*„)). 

Ac  -  C1  ((=,■(»).(  X  ln*n)  *  L2((-r,0>,Ln>j)  . 

Then  Ac  C  A  and  for  \  =  (AQ , A^ • • • , Ay, A)  =  (AQ(* , X) ,A1(* , A) , • • • , 

A  (• ,A) ,A(* , • , A) ) ,  an  element  of  A  or  Ac  define 
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l*L  -  ElAjloo  +  l<  /  |a(*,6) |2de)^|  . 

j=0  J  -r 

Let  L(t;A)  be  the  operator  defined  in  (2.3)  where  the  coefficient  matrices 
are  the  components  of  A.  Let  U^(t,s)  denote  the  solution  operator  and 
x(*,n,<t>,0,A)  the  solution  of  the  FDE  initial  value  problem  (2.1),  (2.2)  with 
f  =  0  and  L(t)  =  L(t;A).  Let  A(t;A):W  C  Z  -*■  Z  denote  the  operator  de¬ 
fined  in  (2.8)  with  L(t)  =  L(t;A)  and  U.  ,(t,s)  denote  the  solution 

A  j  N 

operator  for  the  approximating  initial  value  problem  in  Z^: 

zN(t)  =  AN(t;A)zN(t), 

Vs)  -  ZN 

with 

A„(t;A>  =  PjjAftjA). 

LEMMA  3.1.  Let  A  €  A  be  fixed.  Then  there  exists  a  sequence 
{A^}  C  Ac  such  that  I  Ak  “  ^  k  Moreover,  for  Zq  =  (4>(0),<J>)  €.  W 

we  have 

I UA  (t*s)z0"UA(t’s)zoIz  -  K|Xk"  I z0 1 

k 

where  K  =  K(A)  is  a  constant  independent  of  k  and  Zq. 

Pkoo .  The  existence  of  the  sequence  {A^}  C  Ac  with  A^  -*■  A  is  a 
consequence  of  the  fact  that  A^  is  a  dense  subset  of  A.  Next,  we  let 

(3.7)  z(t)  =  U^(t,s)zQ  =  (x(t;<H0),4>,0,A)  ,xt(<})(0)  ,4>,0,A)) , 

(3.8)  zk(t)  =  U^(t,s)z0  =  (x(t;<HO),<l>,0,Ak),xt(())(0),(f),0,Ak)), 

where  the  extreme  right  hand  equalities  in  (3.7)  and  (3.8)  follow  from  the 
equivalence  established  in  section  2.  We  note  that  for  each  t  ^  s,  z(t),zk(t)  €  W 
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By  [18,  Theorem  3.15],  for  Xk  €.  and  Zq  £  W  we  can  write 

t 

zk(t)  —  Zq  ■+■  £  A(oj X^) 2^(0) do  , 


and  by  [10,  p.  23]  for  X  6  A  and  Zq  £  W 


l 

z(t)  ~  zq  +  f  A(o; X)z(o)do. 


Let  A(t)  =  z(t)  -  z^(t)  •  Then  A(t)  «W  and 


(3.9)  A(t)  =  J [A(o;X)z(o)  -  A(o; Xk)zk(o) ]do 


L 

■  f  jA(o;X) [z(o)  -  zk(o)] +  [A(o; X)  -  A(o;Xk) ]zk(o) jdo 
s  '  ' 

t 

=  J  [A(o;X)A(o) +  6Ak(o)zk(o) ]do. 


where 


6Ak(o)  =  A(0;X)  -A(o;Xk). 


If  we  then  apply  [2,  Lemma  2.1]  (3.9)  implies 


I A ( t ) | 2  =  | A(s) I p  +  2  f  |<A(o;X)A(o),A(o)>f 


+  <6Ak(o)zk(o) ,A(o)>p jdo. 


But  A(s)  =  0  and  (c.f.  Lemma  2.1) 


<A(o;X)A(o),A(o)>p  £  u»(X)|A(o)|  . 


Therefore,  an  application  of  the  Gronwall  inequality  yields 


!  A(t)  1 2  <  2  Jjw(X)|A(0)|2  +  4|6Ak(o)zk(o)|2  +  4|A(o)|2jdo 

-t  j  t  7 

-  J  (2u(X)  +  l)|A(o)rdo  +  J  |5Ak(a)zk(o)rdo, 
s  s 

l/t|6Ak(o)zk(o)l2da  e(2o)(A)  +  l)(t-s)  _ 
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Consider  the  integrand  in  the  last  expression  above: 

r)|2  . 


|6Au(a)2k(a) I' 


=  |([L(a;A)  -  L(a;Ak)]xa((})(0),^,0,Ak),0)  I 

=  I  Y!CA.<a;A)  -  Ai(a;Ak)]x(a-Ti;0(O),<)),O>Ak) 

1=0  1  2 
+  J  [A(a,0;A)  -  A(a,0;Ak)]x(a+e;ij)(0),()>*0,Ak)d«!  • 


Recalling  that  for  Ak  S  Ac  and  zQ 


e  w 


|Ux(t,s)z0i  5  MU0le 

k 


w(Ak)(t-s) 


by  (3.8)  we  have 

)x(T;4>(0),<M,Ak)l  5  Mlz0le 


w(Ak)i 


s-r  <  t  ^  T» 


and  hence 


w(Ak)T 


where  K  is 
for  is  €,  ts>f  J 


|x(*;4>(0) ,4>»0*Ak) |  £  *\*0\*  ~  -  K'z01  ’ 

independent  of  k  for  all  k  sufficiently  large. 


|6Ak(o)zk(o)|2  i  ^U0i  ]£IV,}X)' 

0  1\2 
+  (  J  (a(*,0;A)  -  A(-,e;Ak))de|ot)  j 


^K|z0|  |  Ai<- ;A>  -Ai(-; 


At)  I 


k'  loo 


2 JQ\i I 


+  r 


J  |a(*,0;A)  -  A(*,e;Ak)|  ae) 
-r 


•1)! 
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/  i>  ,  \h  \2 

+  I  (f  |A(*,0;A)-A(*,0;Ak)rde)  Lj 

=  ^2|zoi2  • 

Thus,  we  have 

|uA(t,s)z0  -  UA  (t,s)z0|2  =  ]z(t)  -  z,(t) |2  =  | A(t) | 2 
k 

<  j  /|Mk(o)zk(a)|2da{e<2^  +  1><t-s) 

<  [  /K2|^oi2|A-Aki2d0|e2(ai(X)  +  1)(t~S) 

<  (T-s)K2|z0|2|X-Xk|2  e2(w(A)  +  l)(T-s)t 

which  implies 


luA(t,s)z0- UA  (t,s)z0|  £  K  I  z0  i  lA"Xkloo  » 

k 

where 

K  =  K(T-s)^e^A^  +  1^T-S^  . 


Let  X  =  (aq>Aa>  •  •  •  ,Ay  ,A)  eA  be  given  and  consider  the  FDE  initial  value 


x(t)  =  L(t;X)xt> 
(x(s),xg)  =  zQ  =  (n»4>)  • 


problem 
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Using  the  fact  that  W  is  dense  in  Z  and  A is  dense  in  A,  for 
e  >  0  given,  we  can  choose  z^  €  W  such  that 

| z  -51  <  -  e"u,(A)(T_s) 

1  0  01  M 


and  X  €  A  such  that 
c 


|x-xl 


K(X)(|  e~w(X)(T  s)+  |zQ|)  ' 


Since 


zo  i  lzo“zoi 


+  l*ol  ♦ 


we  have 


ML  < 


K(X) | : 


Now, 


(3.10)  |[PH0A(t,8)-UX,N(t,s)PNlz0|  <  |PNUx(t,s)[z0-i0]| 

+  |[PNUA(t,s)-PNUx(t,s)]Zol  +  |[PNUx(t,s)-Ux>N(t,s)PN]S0| 

+  lUA,N(t’s)PN[i0-z0]l  1  MeW(A)(T-s)|z-Z0|  +  K(X)|2Q|  |X-X|c 

+  I  [pNUx(t,s)  -Uj[jN(t,s)PN]z0l  +  M  |z-i0|eal(A)(T"s) 

<  |  +  K(X)  |  zQ| 

+  |[PNU~(t,s)-Ux>N(t,s)PN]z0l+  Meu3(A)(T_s)  £  e-a,(A)(T-s) 

-  2c  +  |[PNUx(t,s)-U^N(t,s)PN]Z0|  +  EeW)-“W)d-s) 

-*■  0  as  e  -*•  0, 
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where  the  second  term  in  the  final  expression  above  tends  toward  zero 
uniformly  in  t  for  t  e:  [s,T]  by  Theorem  3.1  and  the  coefficient  of  the 
third  term  tends  toward  1  as  f.  -+  0  as  a  consequence  of  the  continuity  of 
w  with  respect  to  A. 

To  summarize,  (3.10)  above  reveals  that  in  order  to  obtain  an  approxi¬ 
mate  solution  to  (2.1),  (2.2)  corresponding  to  A  €  A  it  suffices  to  apply 
our  approximation  schemes  to  an  approximating  FDE  initial  value  problem. 

That  is,  apply  them  to  (2.1),  (2.2)  corresponding  to  \  6  a  smooth 
approximation  to  X.  However,  when  actually  implemented  in  practice,  the 
approximation  schemes  which  we  have  developed  rely  upon  the  application  of 
standard  discrete  numerical  methods  for  ordinary  differential  equations  to 
the  initial  value  problem 

(3.11)  ( t)  =  Ajj(t;  A)z^(t) 

(3.12)  zN(s)  -  PNz0  . 

If  we  replace  A  (t;A)  with  A„(t;A)  and  if  the  time  step  in  the  ordinary 
N  N 

differential  equation  integration  is  chosen  sufficiently  small,  the  resulting 
numerical  solution  would  be  indistinguishable  from  the  one  obtained  by  simply 
integrating  (3.11),  (3.12)  as  it  stands.  Thus,  although  the  convergence 
result  stated  in  Theorem  3.1  applies  only  to  FDE  initial  value  problems  with 
smooth  coefficients,  in  practice  our  approximation  schemes  are  applicable  to 
FDE  with  right  hand  sides  in  L  as  well. 
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4.  Spline  Approximations  and  Their  Implementation 

In  this  section  we  outline  the  ideas  involved  in  realizing  the  approxi¬ 
mation  schemes  discussed  in  section  1  in  a  manner  appropriate  for  computer 
implementation.  The  formulation  employed  in  [8]  for  schemes  developed  for 
autonomous  equations  can  be  modified  so  as  to  be  applicable  to  schemes  for 
nonautonomous  equations  as  well.  Therefore,  with  the  exception  of  the  modi¬ 
fications  required  by  the  time  dependence  of  the  operator  A(t),  the  results 
below  are  a  summary  of  results  found  in  [8].  We  conclude  the  section  with 
the  construction  of  a  particular  realization  using  spline  functions  which 

satisfy  the  hypotheses  and  conditions  of  Theorem  3.1. 

Let  {Z„,P..,A„(t)  }  be  an  approximation  scheme  for  the  FDE  initial  value 
N  N  N 

problem  (2.1),  (2.2)  which  satisfies  the  hypotheses  and  conditions  of  Theorem 

3.1.  Assume  dim  ZN  =  ^  N  »  1,2,***  .  We  recall  condition  (1)  of  Theorem 

N  N 

3.1,  which  states  that  Z„,  C  W,  and  fix  a  basis  for  Z  ,  6.  =  (8.(0), 3.) 

N  N  J  J  3 

j  =  1,2,  •  •  • ,  kN,  with  gj  €H1(-r,0).  Let  8N  denote  the  n  >.  matrix 
function  defined  on  (-r,0)  by 


6N  =  (6”,S?,-*-,8j  ), 

1  z  kN 


and  let 


3N  =  (BN(0),BN). 


For  any  z_.  €  Z..,  we  can  write 
N  N 


ZN  =  =  (BN(0)on,6Non), 


'N 


where  CR  is  the  coordinate  vector  representation  of  z^  with  respect 

-N  k 

to  the  basis  {8. } .  . . 

J  1  =  1 


T-T" 
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Central  to  our  approximation  schemes  is  the  finite  dimensional  approxi¬ 
mating  ordinary  differential  equation  initial  value  problem  in  given  by 


(4.1) 


Vt}  =  AN(t)zN(t)  +  PN(f(t),0), 


(4.2) 


ZN(S)  = 


g 

If  we  let  w^(t),  F^j(t),  and  denote  coordinate  vector  representation  of 

zN(t) ,  PN(f(t),0),  and  PN(n»4>)  respectively,  and  if  we  let  A^(t)  denote 
the  matrix  representation  for  the  operator  A^(t),  all  with  respect  to  the 


basis  the  initial  value  problem  (4.1),  (4.2)  reduces  to 


w^(t)  -  Ajj(t)w^(t)  +  F^(t) 


WN(S)  =  V 


which  can  then  be  solved  via  standard  numerical  methods  for  the  integration 
of  ordinary  differential  equations.  However,  in  order  to  do  this,  we  must 
first  compute 


(l)  PN(n,4>)  €  ZN  for  (n,<J>)  e  z 


(2)  A^t),  t> 


'ith  respect  to  the  basis  We  begin  with  (1). 


Since  P^  is  the  orthogonal  projection  Z^  Z^,  the  orthogonality 


relationship  in  Z 


(4.3) 


{PN(n,<}))  -  (n , 4>) )  |  zN 


uniquely  determines  P^(n»4>)>  and  therefore  w^  as  well.  Expression  (4.3) 
is  equivalent  to 


<8N,^wJ  -  (n ,<b) >p  =  o. 
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which  Implies 


<BN,A;>p  =  <BN,(n,4>)>p  * 


or 


(A.  4) 


V*  -  ’ 


where 


.  <f,f>  =  &n(0)t6n(0)  +  / eN(0)V(e)p(e)de  , 


Qn  p 


and 


Vn,«  -  4",<n,«»0  -  *  /,  6S<»I«e>^9>'le 


Therefore, 

»» ■  %l  V"-',>  ■ 

The  calculations  above  provide  a  Keans  by  "bleb  FN(t>  can  bc  COmra 
well.  Indeed,  (4.4)  implies 

VN(t)  =  hN(f(t)’°)’ 


but 


hM(f(t),0)  =  6N(0)Tf(t), 


and  hence. 


F^(t)  =  QN1BN(0)Tf(t)' 


We  next 


address  (2)*  «t  ♦,  E  Z»  SO,>P°“  “« 


S) 

«  r  is  such 


that 


Furthermore,  for  each  t  -1  s,  let  Yjg(t)  he  such  that 


V'>*«  ■ 

It,  of  course,  then  follows  that 

(4.5)  YN(t)  =  Vt)CV 

Since  ^N(t)$N  =  P/(t)$N  -  PN(L(t)<(>N,D0N),  Y„(t)  is  the  coordinate  vector 

representation  for  (L(t)c}>jj,D<}>^)  .  Therefore,  by  (4.4) 

QNYN<t)  =  VL(t)VDV  =  hN(L(t)BNaN,(DBN)aN)  =  Vt)(V 


where 


iyt)  =  hN(L(t)BN,(DBN)) 

-  eN(o)T(L(t)3N>  +  f  BN(0)T(D3N)(0)P(e)de  . 
-r 


Thus, 

QNYN(t)  =  HN(t)V 

or 

V°  =  QN\(t)aN* 

which  by  (4.5)  implies 

Vc)  ■  ^V0- 

Since  we  have  assumed  that  {Z^,P^,A^(t) }  satisfies  the  hypotheses  and 
conditions  of  Theorem  3.1,  it  follows  that 
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lim  6  w  (t)  =  (x(t;r)»<f>,f)  ,xt(n,(J),f)) , 

N  -*■ 00  ' 


lim  SN(0)wN(t)  =  x(t;n,$,f)  , 


uniformly  in  t  for  t  S  [s,T]  and  uniformly  in  f  for  f  in  bounded 
subsets  of  L^Cs.T). 

We  conclude  this  section  with  the  description  of  a  particular  approxi¬ 
mation  scheme  which  is  included  in  the  framework  constructed  in  section  3. 
While  the  scheme  we  develop  is  closely  related  to  the  spline  based  schemes 
discussed  in  [6]  and  [8],  slight  differences  in  formulation  necessitate  the 
presentation  of  the  scheme's  development  in  detail. 

N  N  Nr 

Consider  the  partition  of  [ — r , 0 ]  given  by  { t ^ ^ with  t^  =  - j  — 
j  =  0 , 1, 2 , • • • ,N,  and  define 

Z1  =  {(r),4>)  =  <{>(0),  <J>  a  first  order  spline  function  with 

N 

knots  at  t^ ,  j  =  0, 1,2, • • • ,N}. 

The  set  is  a  finite  dimensional  subspace  of  Z  with  dim  Z^  =  n(N+l). 

A  basis  for  Z^  may  be  constructed  as  follows: 

For  {tj}j_g  as  above,  and  each  j  =  let  e^  (•)  :  [-r,0]  -*•  R 

denote  the  "hat"  functions  defined  by 


e£(0) 


f  (0-  tj) 


t"  <  e  <  o 


otherwise 


ejN(9)  = 


r  (0  W 


H  /fl  N 
r  (0-tj+l> 


■J  *  •  *  'U 


N  N 

‘j+i  -  0  -  tj  *  J  = 1'2-' 


,N-1  , 


otherwise 
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and  let 


e>  - 


-r  <  6  <  t 


N 

N-l 


otherwise 


6y+i)k  =  (ej  (®)vk»ejvk)  J  ~  0,1,2, ••• ,N  ,  k  =  1,2, , 

T  ti 

with  vk  =  (0, ••• ,0,1,0, ,0)  €  R  where  the  1  appears  in  the  kth  position 

It  can  easily  be  verfied  that  is  a  basis  for  Zjjj.  Let  P*  be 

the  orthogonal  projection  Zp  -*•  Zjjj  and  define  A^(t)  =  pjjjA(t). 


THEOREM  4.1.  The  approximation  scheme  {Z.^,P^,A^(t) }  satisfies  all 
-  N  N  N 

of  the  conditions  in  the  statment  of  Theorem  3.1. 


P>ioo 6 .  Since  conditions  (1)  and  (2)  of  Theorem  3.1  are  trivially 
satisfied  by  {Z^,P^, A^(t) } ,  we  need  only  to  argue  that  it  satisfies  condi¬ 
tion  (3)  as  well.  That  {Z^,P^, A^(t) }  satisfies  condition  (3a)  is  established 
in  the  proof  of  Theorem  4.1  in  [8].  Therefore  we  only  discuss  condition 
(3b)  here. 

Let  $  =  (<j>(0),<t>)  W  with  <p  6H2(-r,0)  and  let  =  P^$  =  (<J>N(0)  ,<J>N)  • 
Theorem  2.5  of  [22]  (see  also  Theorem  21  of  [23])  implies 


(4.6) 
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and 


k-1 


(4.7) 


k-1 


/  i  -h  <§>4  /  i»2*i 


I  2 

where  <j>^  denotes  the  interpolatory  spline  for  <p  €.  H  (~r,0)  with  knots 

at  From  (4.6)  and  (4.7)  we  find 

(4.8)  |D(4>  -  4>*)  1 2  <  ~  |D2<0|2  , 

1  ,r.  2 t_2 ,  | 


N'2  -  *T  (N}  'D  M’ 1 2  * 
n 


(4.9)  |<f>  -  <t> 

Making  use  of  the  norm  equivalence  relation 


lz  1  1*1,  i  ^ 


together  with  the  minimality  properties  of  the  orthogonal  projection 

P2:Z^  Z*  we  find 
N  p  N 


(4.10) 


v-*i2  < 

i  !®n-$Id  i  i^-»i2,0 

<  (£)z  )d2*|  2  , 


where 


t>l  =  (♦£«»,♦£>  =  W(0),4)J)  €  zj 


We  next  use  the  Schmidt  inequality  [22]  to  estimate  jD((J>^  -  cp^.)  |  Since 


N 


J>^,  <t>^  are  linear  on  each  sub-interval  [tj,tj_^]  we  have 


(4.11)  |D4N-o.J)|2  = 


J\v%-*1)\2  -  £  / 


1  fc-c(f)2/’V^|2 

j  =  l  *^N  N 

tN  1  tN  x 

J  '  KT  XI 


iK(f)2  /|V*  I2  +  Vl>2£  / 

J  =  1  N 


|D2^!2 


,N. 2 | .  , | 2  J  k  .rv2._2.i2 

=  k  (-)  I <f>N -  4> |  2  +  -j  (^)  |D  <t>|2 

IT 


.N.2  v  .r.  4 .2  | 2  k  .rv2.  2  ,2 
i  <  (-)  -j  <jj)  |D  4>l2  +  -£  (p  |D  <H2 

TT  TT 


-  (v+1)  ^  (f)2  |D2^2  > 

TT 


where  we  have  used  (4.7)  and  (4.10)  In  making  the  estimates  above.  Therefore 
by  (4.9)  and  (4.11)  we  find 

(4.12)  |d(*h-*)|2  <  !d(4>n-4>J)!2  +  !d(4>J -  4>)  t2 


<  ^n<(L)l^\2  +  j_^l^\2 

TT  TT 


<  K(N)  ID  4)1, 
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where 

K(N)  =  0(|)- 

Noting  that  K(N)  £  K,  K  independent  of  N  and  <j>  6  H2(-r,0),  we  have 

|d(4>n  -  4>)  12  0  N  -*■  °°, 

and 

1d<4>n-  4>)  |2  <  k|d24>|2, 

which  establishes  the  second  part  of  condition  (3b).  To  see  that  for  each  t, 
L(t)cf>N  ■+  L(t)(f>,  it  can  be  argued  (cf.  [8],  Theorem  4.1)  that  (4.12)  im¬ 
plies 

|4>n(9)  -  cf>(0)  |  £  0(i), 

as  N  -*■  00 ,  uniformly  in  0  for  0  S[-r,0].  Therefore,  for  each  t  £  s, 
|L(t)4>N  -  L  ( t )  <p  j  £  0Cj~)  as  N  -*■  °°  and  condition  (3b)  has  been  established. 

For  the  Z*  basis  defined  earlier,  the  matrices  QN  and 

H^t)  take  on  a  particularly  simple  form.  Indeed,  for  the  case  v  =  1 
(and  therefore  p(0)  =  1)  we  have 


(t)  =  f  A(t ,  0)e^(0)d0 

-r  J 


for  N  =  2,3, •••  where  0  denotes  the  Kronecker  product  and  I  is  the 
n  x  n  identity  matrix. 

As  in  the  case  of  the  approximation  schemes  developed  for  autonomous 
equations  in  [8], modifications  of  the  results  presented  above  can  be  used 
to  verify  that  approximation  schemes  employing  higher  order  spline  functions 
and  satisfying  the  conditions  in  the  statement  of  Theorem  3.1  can  be  construe 
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5.  Approximation  Schemes  Incorporating  Time  Discretization 

In  this  section  we  briefly  outline  and  discuss  two  alternative  approxi¬ 
mation  scheme  formulations  based  upon  the  approximation  framework  and  state 
discretization  techniques  developed  in  section  3.  However,  the  two  schemes 
which  we  are  going  to  introduce  differ  from  the  strictly  semi-discrete  schemes 
discussed  previously  in  that  these  schemes  incorporate  various  degrees  of 
time  discretization  together  with  the  state  discretization.  In  particular, 
the  first  alternative  allows  for  the  discretization  of  the  time  dependence 
of  the  operator  L(t).  This  capability  is  extremely  desirable  when  the  time 
dependent  matrix  coefficients  on  the  right  hand  side  of  the  FDE  are  computa¬ 
tionally  expensive  to  evaluate.  On  the  other  hand,  the  second  scheme  allows 
for  the  complete  discretization  of  all  time  dependence  appearing  in  the  equa¬ 
tion.  These  schemes  result  in  difference  equation  approximations  and  are  in 
the  same  spirit  as  the  approximation  techniques  discussed  in  [18]  and  [20]. 

We  discuss  each  alternative  separately  and  In  turn. 

Let  {ZN,PN, A^t) }  be  an  approximation  scheme  for  the  initial  value 

/\ 

problem  (2.1),  (2.2)  and  for  each  t  _>  s  define  the  operators  A^ft^Z^  -*•  Z^ 
by 

(5.1)  VtJ-Vkf).  (k-l)g<t<k 

r  r 

for  k  =  i,  i+1,*’*,  where  i  is  that  integer  for  which  (i-1) —  <  s  £  i  — . 

We  consider  the  ordinary  differential  equation  in  Z^  given  by 

(5.2)  zN(t)  =  AN(t)zN(t). 

Since  A^(t)  *-s  piecewise  constant,  the  evolution  operator  corresponding  to 

(5.2)  is  of  the  following  form: 


Zl'l 
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If  (i-D  f  <  s  1  c  1  1  §  . 


UN(t,s)  =  exp  |(t-s)A^(i  J)j. 

If  (i-1)  l  <  s  <i  §<  •••  <j-l  f  <  t  <  j  §, 

UN(t,s)  =  expj(t-  (j-1)  f>VJN)]eXp  §Sj((j~1)  N5]  * 

. eXP[^((i+1)i)]  6XP  [(±f  '  s)AN(if>]- 


It  then  follows  that 

-  .  <  U,(t-(^1)i)  W(iN~s)  “<f>  (t-s)a) 

|UN(t,s) (  <  e  e  | I  e  =  e 

N  P  k=i+l 


and  therefore  for  t  <  T 


|UN(t,s)i  <  Mew(T“s), 


where 


M  =  /  V  . 


THEOREM  5.1.  Consider  the  FDE  initial  value  problem  (2.1),  (2.2)  under 

the  additional  assumption  stated  in  Theorem  3.1.  Suppose  further  that 

{Z„,P„T,A  ft) }  is  an  approximation  scheme  satisfying  conditions  (1),  (2),  (3) 

N  N  N 

of  Theorem  3.1.  Then  if  (Z^,PN,A^(t) }  is  an  approximation  scheme  with  A^(t) 
defined  as  in  (5.1),  we  have 

lim  |[P  U(t,B)-U  (t,s)P  ]zQ|  =  0  , 

N  -*■  °° 

for  each  Zg  £  Z  uniformly  in  t  for  t£[s,T]< 
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PMO±.  As  in  the  proof  of  Theorem  3.1,  It  can  be  argued  that  for 


z06D 


|UN(t>s)PNZ0~PNU(t’s)zol  =  I  /UN(t,T)[XN(T)PN-PNA(T)]U(T)s)z0dx| 

s 

T 

^  M  W(T— S  )  <*  I  r  A  ..  t>  A  / \  1.  /  I 


Me  J"  |  [\(t)pn- PnA(t)]z(t,s)  | dt  . 


The  desired  result  will  clearly  follow  if  we  can  demonstrate  that  for  each 
xG[s,T]  the  integrand  in  the  last  expression  in  (5.3)  above  tends  to  zero 
as  N  with  the  convergence  being  dominated.  Fortunately, 

however,  this  can  be  argued  in  precisely  the  same  manner  in  which  it  was 
done  in  the  proof  of  Theorem  3.1  with  one  minor  exception.  We  must  show 
that  for  (j-l)~  <  T  <_  j|,  z(t,s)  =  (x(t),xt>  and  Pnz(t,s)  E  (xN(t),x^), 

L(j|)x^  +  L(t)xt  , 

with  the  convergence  being  dominated.  Using  the  estimates  computed  in 
proving  the  analogous  claim  in  Theorem  3.1,  we  have 


|h(j^)x”-  L(t)xt|<  |L(j^)  (x^- xt)  |  +  I  (L(j§)  "  L(t))xtI 


EJAjL  +  r*U  /  ]A(-,e>!2d0)4|co) 

•  (lxN- XL  +  ri|D(x”- xT)  |2) 

/v  0  \ 

+  |a(jJ)-a(t)|  +  f  |A(jJ,e>-A(T,e)|deJ|xT 

i  M1|ixN_xL +  riiD(x^-xT)i2| 

/  v  0  \ 

+  M2^  £  |A(jJ)  -  A(t)  I  +  f  |A(j^,0)  -  A(T,0)|dej 


-*•  0  as  N  -> 
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where  we  have  used  conditions  (3a)  and  (3b)  to  conclude  that  the  first 
term  above  tends  to  zero  and  the  continuity  of  the  right  hand  side  of 
the  FDE  to  conclude  that  the  second  term  tends  to  zero  as  well.  Recalling 
that  by  condition  (3b)  we  have 

|d(x”-xt)|2  <  k|d2xt|2  <  k|d2x|2> 

it  then  follows  that  the  convergence  in  (5.4)  is  dominated. 

Finally,  as  was  the  case  in  Theorem  3.1,  D  dense  in  Z  and  the  opera¬ 
tors  [UN(t,s)P^  -  P^U(t,s) ]  uniformly  bounded  are  sufficient  to  guarantee 
convergence  on  all  of  Z  and  the  theorem  is  proven. 

We  next  turn  our  attention  to  the  second  alternative  which  involves  the 

total  discretization  of  the  time  dependency  in  the  equation.  Let 

{ZN,pN,AN(t) }  be  an  approximation  scheme  for  the  initial  value  problem  (2.1), 

(2.2)  satisfying  the  hypotheses  and  conditions  of  Theorem  3.1.  Suppose 

further  that  A^(t)  and  i  ate  as  they  were  defined  in  (5.1).  Then  for 

each  N  =  1,2,***,  and  each  t  >  s  define  the  operator  V„(t,s):Z„  -*■  Z„ 

—  N  N  N 

by: 

VN(t,s)  =  jl-  (t-s)AN(t))~1 

i  r 

(i-l)  ~  <  s  <  t  <  i  f 

and 

. 

. (i-sv<i«>s>r1(i-<iR  -‘JV')"1 

if 
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In  light  of  the  arguments  regarding  the  maximal  dissipativeness  of  the 
operators  A^(t)  -  wl  given  in  the  proof  of  Theorem  3.1,  we  have  that  for 
N  sufficiently  large,  the  operator  inverses  required  in  the  definition  of 


VN(t,s)  above  exist.  Furthermore,  using  (3.1)  it  is  easily  verified  that 


TT  -jf 

1 _ /" 


(f, 


.  N 


,  .  /H\  N  N 

k-i  (“)  "  w  as  N 


and  hence,  for  s  t  £  T,  V^(t,s)  is  uniformly  bounded.  In  fact,  it  can 
be  shown  that  for  all  N  sufficiently  large  there  exist  an  to  >  0  such 
that 


(5.5) 


I  Vji  ( t  * s  )  I  1  Me 


io(T-s) 


where 


M  =  /  V. 


THEOREM  5.2.  Let  (Z„,P..,A„(  t)  }  be  an  approximation  scheme  for  the 
-  N  N  N 

initial  value  problem  (2.1),  (2.2)  under  the  additional  smoothness  assumption 
stated  in  Theorem  3.1.  Suppose  further  that  A^(t) }  satisfies 

conditions  (1),  (2),  and  (3)  of  that  theorem.  Then  for  VN(t,s)  as  defined 
above,  we  have 

lim  | [P  U(t,s)  -  VN(t,s)PN]z0|  =  0 

N  ->  co 

for  each  zQ  €  Z  uniformly  in  t  for  t  €  [s,T]. 


The  proof  of  Theorem  5.2  can  be  argued  in  a  manner  similar  to  that 
used  by  Yosida  in  the  proof  of  Theorem  XIV. 2  in  [25].  We  shall  outline  the 
essential  ideas.  Using  the  definition  of  V^(t,s)  and  the  resolvent 
identity  (i.e.  for  A,  B  linear  operators  and  A  6.  p(A)n  p(5)  we  have 
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R(A;A)  -  R(X;8)  =  R(X;A) (A-8)R(X;8)  it  is  not  difficult  to  show  that 
for  N  sufficiently  large 

<5.6)  fj  yt,.)  - 

This  in  turn  implies  that  for  D 

(5-7)  IF  tVN(t*T)PNPNU(T*s)zO] 

*  |l?  VN(t*X)PN  |  PNU(T>S)Z0  +  VN(t>T)PNPN  I?  U(T’S)Z0 

-  Vt.T)|v<T)-VT)  (x-(i  f-T)VT))"lpN|u(T’s)zo- 

Integrating  both  sides  of  (5.7)  from  s  to  t  we  find 

[PNU(t,s)-VN(t,s)PN]z0 

t 

=  f  VN(t,T) |PNA(t)  "  (I_  (i  N"  T)AN(T))'lAN(T)PN  u(T.s)z0dT. 

S 

Using  the  bound  given  in  (5.5)  we  have  for  Zq  €  D  and  z(t,s)  *  U(t,s)z ^ 

|[PNU(t,s)-VN(t,s)PN]z0| 

t 

</lVt.T)|  l(r-  (i  r^VT))'1!  ||AN(T)PN-FNA(T)|z(T,S)|dT 

s 

t 

+  f  |VN(t,x)|  (  [(I-  (i  l|pNA(T)z(T,s)|dT 

s 

T 

<  M2e“(T-s)(l-w  J)"1  f  |[AN(-t)PN-PNA(T)]z(T,s)|di 
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+  M ew(T_s)y*  I  |(l-  (iJ-i)AN(T))-1- I  P^OOzOr.sHdT 

s 

T 

*1  f  I  |Vt)pn'  p/(t)  z(t’s)  ldT 

s 

T 

+  Y2  f\  [  (I  -  (i§-T)AN(T))_1PN-PNl|A(T)Z(T,s)|dT 


T1  +  V 


We  have  already  argued  (cf.  proof  of  Theorem  5.1)  that  ■+  0  as  N  -*■  oc, 
uniformly  in  t  for  t  S[s,T].  Consider  the  term  T^.  Since 

(I_  =  Pl,0  N  ~  T)\(T))  * 

where  P^  g(z)  denotes  the  (1,0)  entry  in  the  Pade  table  of  rational 
function  approximations  to  the  exponential  we  can  apply  [20,  Theroem  10.3] 
and  thus  conclude  that  for  each  T  £[s,T] 


I_  (i  1PN- PNI  A(t)z(t,s)  |  0 


as  N  -*■  °°.  Furthermore,  since  Zg  6  D  and  z(x,s)  =  ,x^)  we  have 


I  -  (i  t)An(t))_1Pn-  PNI  A(t)z(t,s)  | 


I-  (ij  -t)An(t))"1Pn-PniJ|2|A(t)z(t,s)|- 


1(m(1-w^)  1+l)2|  (L(t)xt,Dxt)|2  £  |M(1-W^)  1+i)2  |l(t)xt|2  +  |Dxt 
<  (Md-u.jr^l)2  (J:iAj|<x)+r^|(/(A(.,e))2dG)^L)|xL  2 


+  Dx 


L2(s-r,T) 
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for  N  sufficiently  large.  Therefore  ®  as  N  -*•  ai  uniformly  in  t 

for  t  e [s,T]  and  hence 

|[PNU(t,s)  -VN(t,s)PN]z0|  -  0  as  N  -  »  , 

uniformly  in  t  for  t  e  [s,T]  for  each  £  D.  Since  D  is  a  dense 
subset  of  Z  and  the  operators  [P^U(t,s)  -  V^(t,s)P^]  are  uniformly  bounded 
we  once  again  can  extend  the  convergence  to  all  of  Z  and  the  theorem  is 
proven. 

If  we  define 

zN<t,s;n,0,f)  =  UN(t,s)PN(n,<)>)  +  J  UN(t,T)PN(f(T),0)dT, 

s 

and 

f 

wN(t,s;n,4>.f)  =  VN(t,s)PN(n,<J:)  +  JvN(t,T)PN(f(T),0)dT, 

s 

then  all  of  the  statements  in  Theorem  3.2  regarding  the  convergence  of  the 

approximations  to  the  solution  of  the  nonhomogeneous  problem  remain  valid 

with  z„  replaced  by  either  z„  or  w„T. 

N  N  N 

R-QjrwxAk.  The  operators  U^(t,s)  and  Ujj(t,s)  defined  previously 
must  be  computed  indirectly.  That  is,  they  are  computed  as  the  numerical 
solution  to  the  ordinary  differential  equation  initial  value  problems  in 
ZN  given  by 

UN(t>s)  =  AN(t)UN(t,s) 

UN(s,s)  =  I 


and 
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UN(t.s)  =  AN(t)UN(t,s), 

UN(s,s)  =  I  , 

respectively.  The  operator  V^(t,s)  on  the  other  hand,  can  be  computed 

directly.  Although  in  practice  we  never  actually  compute  the  operators 

U„,  U„  or  V„,  the  above  observation  reveals  that  when  one  uses  the  fully 
N  N  N 

discrete  scheme,  there  is  no  further  approximation  required  beyond  the 
approximation  of  A(t)  by  /^(t)  and  zQ  by  PNzQ.  However,  the  time 
discretization  employed  in  the  definition  of  VN(t,s)  is  essentially  a 
backwards  Euler  scheme  which  is  only  first  order  convergent  (cf.  [20]). 

It  is  unlikely,  therefore,  that  the  discrete  scheme  we  have  presented  would 
perform  as  well  as  the  semi-discrete  schemes  used  in  conjunction  with  say 
a  fourth  order  Runge  Kutta  method  to  integrate  the  approximating  ordinary 
differential  equation.  The  discrete  scheme  should  not  however  be  totally 
discarded  in  that  it  does  represent  a  starting  point  for  the  development 
of  a  discrete  approximation  framework  (similar  to  the  one  developed  for 
autonomous  delay  systems  in  [20])  encompassing  time  discretizations  of 
arbitrarily  high  order. 

6.  Numerical  Results 

In  this  section  we  discuss  a  variety  of  numerical  examples  which  demon¬ 
strate  the  feasibility  of  the  approximation  methods  developed  in  sections 
3  and  4.  All  computations  were  performed  using  a  software  package  written 
for  the  IBM  370  model  158  computer  at  Brown  University.  Since  our  primary 
objective  was  to  test  our  approximation  methods  and  their  convergence  pro¬ 
perties,  factors  such  as  computational  efficiency  and  storage  requirements 
were  given  only  secondary  consideration  in  the  design  of  our  programs.  We 
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note  that  software  packages  developed  to  Implement  the  spline  approximation 
schemes  for  linear  autonomous  functional  differential  equations  discussed  in 
[8]  can  easily  be  modified  so  as  to  be  able  to  generate  approximate  solutions 
to  nonautonomous  equations  via  the  schemes  which  we  have  presented  above. 

In  all  of  the  examples  which  follow,  the  approximation  scheme  employed 
represents  a  realization  of  the  linear  spline  scheme  discussed  in  section 
4.  In  order  to  analyze  the  convergence  properties  of  the  scheme  we  have 
computed  either  exact  solutions  via  the  method  of  steps  or  highly  accurate 
approximate  solutions  using  the  method  of  steps  combined  with  a  fourth  order 
Runge  Kutta  integration  scheme  for  ordinary  differential  equations.  In  the 
latter  case,  we  emphasize  that  although  our  reference  solutions  are  approxi¬ 
mate,  they  were  computed  using  methods  which  are  completely  independent  of 
the  approximation  schemes  which  we  are  testing,  and  hence  should  not  lead 
to  invalid  conclusions.  In  fact,  for  the  examples  for  which  exact  solutions 
were  available,  we  also  computed  approximate  reference  solutions  and  found 
them  to  be  virtually  indistinguishable  from  numerical  tabulations  of  the 
exact  solutions.  The  examples  for  which  the  exact  solution  was  used  as 
the  reference  solution  are  those  for  which  the  appropriate  formulas  have 
been  provided. 

The  entries  6^  in  the  tables  which  follow  represent  the  absolute 
differences  between  the  reference  solution  x  (exact  or  approximate)  and 
our  approximate  solutions  evaluated  at  sample  points  (tk)  along  the  interval 
of  interest.  That  is 

6n  =  |Xj(tk>-  [8N(0)wN(tk)]j|, 

N  N  NT 

for  j  »  1  or  2,  dependinr  on  the  usage  below,  where  w  (t)  =  (w, (t) , • • • ,w„. , (t) ) 

1  N+l 

is  the  solution  vector  of  the  approximating  ordinary  differential  equation 
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N 

computed  using  a  fourth  order  Runge  Kutta  method  and  (3  denotes  the 
matrix  function  on  (-r,0)  defined  in  terms  of  the  linear  spline  basis 
given  in  section  4. 

Example  6.1. 

For  our  first  example  we  consider  the  first  order  homogeneous  equation 
x(t)  =  x(t)  +  tx(t-l), 

on  the  interval  [0,2]  with  constant  initial  data  given  by 

xo(0)  El,  -1  <  0  £  0. 

For  this  example  we  have  computed  the  exact  solution  and  it  is  given  by 

{26*  -  t  -  1  t  £  [0,1] 

2et  +  (t2  -  8)et_1  +  t2  +  2t  +  2  t  £[1,2] 

Upon  inspection  of  the  numer4  'll  results  given  in  Table  6.1,  it  is  easily 
seen  that  convergence  is  second  order. 


Table  6.1 

rt 

x(t) 

** 

m 

6S 

*16 

632 

0 

1.0 

0 

0 

0 

0 

0 

.25 

1.33805 

.0310 

.0108 

.0024 

.0006 

.0002 

.5 

1.7974 

.0655 

.0167 

.0043 

HI 

.0002 

.75 

2.4840 

.1079 

.0275 

.0069 

89 

.0004 

1.0 

3.43656 

.1482 

.0391 

.0097 

.0024 

.0005 

1.25 

4.7773 

.1933 

.0519 

.0145 

.0036 

.0009 

1.5 

6.73324 

.2859 

.0823 

.0205 

.0051 

.0013 

1.75 

9.6190 

.4560 

.1282 

.0325 

.0081 

.0020 

2.0 

13.9050 

.7405 

.2007 

J 

.0516 

.0130 

.0033 
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Example  6.2. 


For  our  next  example  we  consider  the  governing  equation  from  an  optimal 
control  problem  discussed  in  [18] 


x(t)  =  6tx(t-l)  +  u  (t) 


t  €  [0,2] 


xq(6)  =  1 


-1  <  0  <  0 


u*(t)  = 


-3at  -6at+10a 


t  S  [0,1] 


t  €[1,2], 


where  a  =  -23.5/44.8.  The  nonhomogeneous  term  u  appearing  in  the 
equation  above  is  the  point  in  L2(0,2)  where  the  functional 

2 

*<«)  =  Jx2(2)  +  i  f  u2(t)dt, 

Jo 

achieves  its  minimum  subject  to  the  equation  above.  The  exact  solution 
is  given  by 


x(t)  = 


-at3  +  3(l-a)t2  +  lOat  +  1 


(-1.2a)t5  +  4.5t4+  (26a-12)t3 
+  (12-36a)t2  +  at+  (16.2a-  .5) 


t  €  [0,1] 


t  «  [1,2]  . 


The  numerical  results  for  this  example  can  be  found  in  Table  6.2.  Once  again. 


the  convergence  is  essentially  second  order. 
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Table 

6.2 

t 

x(t) 

..  . 

64 

00 

616 

<532 

0 

■ 

0 

0 

0 

0 

.25 

.0238 

.0034 

.0010 

.0002 

.5 

.0002 

.0008 

.0002 

.0000 

.75 

-.14017 

.0089 

.0013 

.0009 

.0002 

1.0 

.85267 

.1353 

.0405 

.0122 

.0039 

1.25 

1.43499 

.1222 

.0095 

.0048 

.0004 

1.5 

.7359 

.1716 

.0082 

.0042 

.0013 

1.75 

-.2029 

.3263 

.0670 

.0125 

.0033 

2.0 

.52455 

.2170 

.0657 

.0186 

.0055 

Example  6.3. 

Here  we  consider  the  scalar  second  order  homogeneous  equation 
x(t)  +  tx(t-l)  +  x(t)  =  0. 

On  the  interval  [0,3]  with  initial  conditions  given  by 


xQ(0)  =  cos (6+1) 


x  (6)  =  -sin (0+1) 


-1  <  0  <  0  . 


Rewriting  the  initial  value  problem  above  as  ^  first  order  system  we  have 

y(t-l) , 


"  0  l" 

’  0  o" 

y(t)  « 

y(t>  + 

I 

1 

,  ° 

0  -t 

V0) 


cos (0+1) 
-sin(8+l) 


-1  <  0  <  0, 


where 


y(t)  =  . 


Our  numerical  results  for  y^  =  x  and  y2  =  x  are  tabulated  in  Tables 
6.3  and  6.4  respectively.  For  this  two  dimensional  example,  the  convergence 
is  observed  to  be  quadratic  for  both  components  of  the  solution. 


Table  6.3 


t 

yx(t) 

0 

.54030 

.5 

.075816 

(H 

2.0 

-.14448 

2.5 

.487092 

3.0 

.86544 

62 

O) 

68 

616 

.00168 

.0003 

.0000 

.0000 

.0046 

.0009 

.0002 

.0001 

.0277 

.0076 

.0019 

.0005 

.0637 

.0163 

.0041 

.0009 

.0885 

.0546 

1 

.0051 

.0018 

.0012 

.0003 

.0868 

.0291 

.0079 

.0022 

Table  6.4 


y2<t) 

62 

-.84147 

.0021 

-.95738 

.0126 

-.62366 

.0247 

.17833 

.0221 

1.07236 

.0198 

1.2591 

.1283 

.01298 

.2552 

.0003 

.0035 

.0060 

.0050 

.0081 
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Examp le  6.4. 

We  again  consider  a  damped  oscillator  with  delayed  damping.  In  this 
example,  however,  we  place  the  time  varying  coefficient  in  front  of  the 
restoring  term  instead  of  the  damping  term  as  was  the  case  in  the  previous 
example.  The  initial  value  problem  is  given  by 


x(t)  +  x(t-l)  +  tx(t)  =0  t  £ [0,3]  , 

x  (0)  =  cos (8+1) 

-i  <  e  <  o, 

xQ(0)  =  -sin(0+l) 

y(t-l)  t  6 [0,3]  , 


0  1 

ro  oi 

y(t)  = 

-t  o_ 

y(t)  +  ' 

1 

- 1 

T—\ 

1 

0J 

where 


y0(6) 


COS (0+1) 
-sin(0+l) 


y(t) 


x(t)- 

x(t)_ 


-i  £  e  £  o  , 


if  written  as  an  equivalent  first  order  system.  Examination  of  the  numerical 
results  contained  in  Tables  6.5  and  6.6  reveals  that  the  quadratic  convergence 
is  unaffected  by  the  placement  of  the  time  varying  coefficient.  Once  again 
we  have  second  order  convergence  in  both  the  x  and  x  components. 


Table 

3.5 

ft 

yl(t) 

64 

68 

616 

0 

.540302 

1 

.5 

.133213 

.0121 

1.0 

-.166604 

.0096 

.0006 

1.5 

-.230397 

.0152 

.0037 

.0009 

2.0 

-.035350 

.0480 

.0098 

.0023 

.0005 

2.5 

.262053 

.0240 

.0023 

.0006 

3.0 

.365969 

.1226 

.0307 

.0076 

.0019 

Table 

6.6 

t 

y2(t) 

6, 

4 

58 

616 

-.841470 

.0021 

.0003 

.0000 

.0000 

-.752190 

.0184 

.0044 

.0012 

.0003 

-.396825 

.0016 

.0004 

.149976 

.0012 

.0004 

2.0 

.577757 

.0884 

.0056 

.0015 

2.5 

.504380 

.1601 

.0384 

.0094 

.0024 

3.0 

-.163431 

_ 

.0899 

.0148 

— 

.0029 

_ . _ 

.0006 

Example  6.5. 

In  this  example  we  consider  a  damped  oscillator  with  sinusoidal 
external  force  and  delayed  restoring  force  having  a  time  varying  coefficient. 
The  initial  value  problem  is  given  by 


x(t)  +  x(t)  +  tx(t-l)  =  sint 
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or 


where 


xQ(0)  =  sin (0+1) 
xq(G)  =  cos (0+1) 


-1  <_  0  <_  0, 


o  r 

"o  o' 

"  o' 

y(t)  + 

y(t-i)  + 

0  -1 

L-t  oj 

sin  t 

y0(9) 


s  in  (0+1)' 
cos(0+i)_ 


-1  <  9  <  0  , 


y(t) 


x(t)' 

x(t) 


As  the  numerical  results  in  Tables  6.7  and  6.8,  substantiate,  second  order 
convergence  in  both  the  x  and  x  components  of  the  solution  is  maintained 
even  with  the  inclusion  of  the  nontrivial  forcing  term.  However,  as  will 
become  apparent  from  the  numerical  results  discussed  in  Example  6.9,  the 
rate  of  convergence  is  sensitive  to  the  smoothness  of  the  function  appearing 
as  the  nonhomogeneous  term  in  the  equation. 


Table 

>.7 

t 

✓"N 

rt 

w 

62 

64 

68 

616 

0 

.841470 

.0021 

.0003 

.0000 

.0000 

.5 

1.067646 

.0006 

.0007 

.0001 

.0000 

1.0 

1.243813 

.0125 

.0027 

.0007 

.0002 

1.5 

1.341552 

.0152 

.0040 

.0003 

2.0 

1.269265 

.0057 

.0006 

.0000 

2.5 

.902306 

.0147 

.0049 

.0003 

3.0 

.132715 

.0425 

,0115 

_ 

.0007 
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Table  6.8 

t 

y2(t) 

*2 

54 

68 

516 

0 

.540302 

.0017 

.0002 

.0000 

.0000 

.5 

.395912 

.0174 

.0034 

.0008 

.0002 

1.0 

.296513 

.0029 

.0008 

1.5 

.063469 

.0033 

.0007 

.0004 

.0001 

2.0 

-.395815 

.0115 

.0047 

.0009 

.0002 

2.5 

-1.111002 

.0238 

.0054 

.0014 

.0003 

3.0 

-1.97221 

.0355 

.0091 

.0021 

.0005 

Example  6.6. 

For  the  second  order  equation 

x(t)  +  e  tx(t-l)  +  x(t)  =  sin t  t  e  [0,3]  , 

with  initial  data 


XqCO  =  cos  (3+1) 


xQ (0)  =  -sin (0+1) 


-1  <  0  <  0  , 


the  numerical  results  given  in  Tables  6.9  and  6.10  exhibit  second  order 


convergence . 
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Table 

j.9 

t 

x(t) 

62 

n 

68 

516 

0 

.540302 

.0017 

.0002 

.0000 

.0000 

.5 

.106974 

.0093 

.0026 

.0002 

1.0 

-.171978 

.0413 

.0108 

.0007 

1.5 

-.134712 

.0671 

.0167 

.0010 

2.0 

.217952 

.0554 

.0123 

.0029 

.0007 

2.5 

.745662 

.0035 

.0033 

.0010 

.0003 

3.0 

1.227099 

.0904 

.0239 

.0061 

.0016 

Table 

3.10 

t 

i(t) 

62 

64 

oo 

<o 

616 

0 

-.841470 

.0021 

.0003 

.0000 

.0000 

.5 

-.791286 

.0217 

.0061 

.0015 

.0004 

1.0 

-.269420 

.0232 

.0050 

.0012 

.0003 

1.5 

.416949 

.0127 

.0039 

.0011 

.0003 

2.0 

.944983 

.0766 

.0206 

.0052 

.0013 

2.5 

1.08968 

.1262 

.0308 

.0077 

.0019 

3.0 

.756316 

.1253 

.0284 

.0069 

.0017 

Example  6.7. 

In  this  example  we  consider  the  same  damped  oscillator  as  the  one  dis¬ 
cussed  in  Example  6.5.  Here,  however,  we  exclude  the  external  forcing  func¬ 
tion  and  provide  discontinuous  initial  data.  The  initial  value  problem  is 
given  by 


x(t)  +  x(t)  +  tx(t-l)  =  0 


t  e  [0,3] 
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x0(e) 


I 


1 


-1 


e  c(-J,o] 


X 


0 


(0)  =  0. 


The  numerical  results  for  this  example  can  be  found  in  Tables  6.11  and 

6.12.  From  the  data  contained  In  these  tables,  it  can  easily  be  inferred 

that  the  rate  of  convergence  is  sensitive  to  the  smoothness  of  the  function 

given  as  initial  data.  Indeed  the  convergence  exhibited  in  both  the  x  and 
* 

x  components  of  the  solution  is  clearly  not  second  order.  It  is  interesting 
to  note  initial  data  that  is  continuous  but  not  can  also  lead  to  sub¬ 

quadratic  convergence.  We  shall  see  this  in  the  next  example. 


Table  6.11 

t 

X(t) 

62 

64 

*8 

616 

0 

-1.000 

.0430 

.0070 

.0000 

.5 

-1.018469 

.0094 

.0324 

.0260 

.0069 

-.994824 

.0482 

.0549 

.0211 

.0013 

-.779472 

.0701 

.0112 

.0137 

.0010 

2.0 

-.344241 

.1063 

.0615 

.0164 

.0004 

2.5 

.304900 

.1421 

.0257 

.0064 

.0041 

3.0 

1.063963 

.1508 

.0018 

.0028 

.0059 
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Table  6.12 

t 

*(t) 

62 

64 

68 

516 

0 

0.0 

0 

0 

0 

.5 

-.106530 

.0766 

.0442 

.0143 

.228874 

.0286 

.0349 

.0114 

.0120 

H 

.641839 

.0407 

.0411 

.0263 

.0033 

2.0 

1.099665 

.0880 

.0161 

.0167 

2.5 

1.463359 

.0484 

.0838 

.0074 

.0034 

3.0 

1.494070 

.0533 

.0387 

.0018 

.0066 

Example  6.8. 

Here  we  consider  the  second  order  equation 

x(t)  +  tx(t-l)  +  x(t)  =0  t  6  [0,3], 

with  continuous  but  not  initial  data 


i  +  e  e  e  [-i,-i ] 


xo(0) 


j1 

ve>  ■  L 


Written  as  a  first  order  system, 


9  fe[-4,  0] 

e  €[-1.-41 

0  e.  (-4.  o] 

the  initial  value  problem  becomes 


"0  1  ' 

0  o  ‘ 

rt 

V-' 

H 

y(t)  + 

0  _ 

o  -t 

y(t-i). 


1 
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with 


y(t) 


x(t)‘ 

i(t)_ 


s  e  [-1,-1] 


e  e (-i,o] 


It  is  evident  from  Tables  6.13  and  6.14  that  second  order  convergence  has  not 
been  achieved.  In  fact  there  does  not  appear  to  be  a  discernible  pattern 
to  the  convergence.  It  can  be  pointed  out  that  the  initial  data  for  this 
example  is  contained  in  ,  N  2,  where  denotes  the  linear  spline 

approximation  spaces  being  used  (see  section  4).  It  is  evident  from  this 
example  that  this  feature  need  not  have  a  positive  effect  upon  the  performance 
of  the  approximation  schemes.  This  is  in  contrast  to  the  behavior  observed 
for  the  averaging  schemes  discussed  in  [8]  when  applied  to  autonomous  equations. 


Table  6.13 

t 

yx(t) 

62 

64 

07 

00 

616 

0 

0.0 

0.0 

0.0 

0.0 

.5 

-.5 

.0361 

.0122 

.0045 

1.0 

-.843932 

.2110 

.0138 

.0149 

.0113 

1.5 

-.736574 

.3184 

.0264 

.0263 

.0202 

2.0 

-.123418 

.2915 

.0013 

.0201 

2.5 

.629872 

.0372 

.0027 

.0108 

.0074 

3.0 

.756091 

.3846 

.0384 

.0001 

.0105 
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Table  6 

14 

t 

y2(t) 

52 

6 

4 

68 

616 

_ 

0 

-1.0 

.0819 

.0122 

■MB 

.5 

-1.0 

.0735 

.0717 

■ 

I 

1.0 

-.289816 

.1923 

.0159 

.0159 

1.5 

.752107 

.0427 

.0296 

.0091 

2.0 

1.57099 

.2871 

.0718 

.0383 

.0180 

2.5 

1.172884 

.6506 

.0309 

.0580 

.040 

3.0 

-.920864 

.6728 

.0716 

.0309 

.0239 

Examp  le  6 . 9 ■ 

For  our  final  example  we  again  consider  the  damped  oscillator  with 
delayed  restoring  force  of  Example  6.5.  Here,  however,  we  include  a  dis¬ 
continuous  external  forcing  function  in  the  equation.  The  initial  value 
problem  is  given  by 


x ( t )  +  x ( t )  +  tx(t-l)  =  10u  ^(t)  t  e  [0,3], 


xQ(0)  =  cos(0+l) 
x0(6)  =  -sin (0+1) 


-i  <  e  £  o, 


where  u  denotes  the  unit  step  at  t  =  .5  defined  by 


u  5(t)  = 


t  <  .5 
.5  <  t. 


Rewritten  as  a  first  order  system,  the  above  problem  becomes 


-56- 


y(t)  = 


0 

1  ■ 

no  o' 

V 

0 

-Ij 

y(t)  + 

-t  Oj 

y(t-i)  + 

10 

u.5(t)  1 


y0(6) 


cos (0+1) 
-sin (0+1) 


-i  <  e  <  o, 


where 


y(t)  = 


x(t  )" 
x(t)_ 


Not  unexpectedly,  we  found  that  the  observed  rate  of  convergence  for  our 
approximation  schemes  was  sensitive  to  the  discontinuity  in  the  forcing 
function.  As  the  numerical  results  in  Tables  6.15  and  6.16  bear  out, 
convergence  in  the  first  component  of  the  solution  appears  to  be  second 
order  throughout,  while  in  the  case  of  the  x  component,  quadratic  con¬ 
vergence  is  obtained  at  some  of  the  sample  points  but  not  at  others.  The 
discontinuity  in  the  nonhomogeneous  term  introduces  a  jump  in  the  second 
derivative  of  the  solution.  This,  in  turn,  will  degrade  the  performance 
of  the  approximation  schemes. 


Table  6.15 

t 

yL(t) 

62 

64 

68 

616 

0 

.540302 

.0017 

.0002 

.0000 

.5 

.191454 

.0182 

.0044 

B 

.0003 

1.0 

.982577 

.1672 

.0454 

.0027 

1.5 

3.330325 

.2486 

.0585 

H 

.0037 

2.0 

6.631893 

.2137 

.0551 

.0133 

.0031 

2.5 

10.105261 

.0876 

.0077 

.0017 

.0001 

3.0 

12.380963 

.2126 

.0654 

.0177 

.0050 

-57- 


Table  6 

16 

t 

y2(t) 

62 

64 

68 

616 

0 

-.841470 

.0021 

.0003 

.0000 

.5 

-.557825 

.0214 

.0043 

.0023 

.0010 

3.386594 

.1169 

.0115 

.0047 

.0012 

E9 

5.823277 

.0522 

.0010 

2.0 

7.133736 

.0543 

.0261 

.0104 

.0036 

2.5 

6.298532 

.2048 

.0628 

.0143 

.0041 

3.0 

2.180855 

.4984 

.1208 

.0318 

.0086 

Although  a  complete  characterization  of  the  convergence  properties  of  the 
approximation  schemes  we  have  developed  would  be  extremely  difficult  if  not 
impossible  to  obtain,  based  on  the  numerical  evidence  which  we  have  presented, 
the  following  conclusion  can  safely  be  drawn.  Our  schemes  appear  to  yield 
quadratic  convergence  for  initial  value  problems  in  which  the  initial  data, 
external  forcing  function,  coefficients  and  therefore  the  solution  are 
smooth.  On  the  other  hand,  it  is  a  relatively  simple  matter  to  break 
the  second  order  convergence  in  some  or  all  of  the  components  of  the 
solution  through  the  introduction  of  irregularity  into  the  solution. 
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